#!/usr/bin/env python
import matplotlib
matplotlib.use('Agg')
import pylab
import re
import argparse
import sys
from scipy.ndimage.filters import gaussian_filter1d
class bandProp():
    def __init__(self,title):
        self.eneRange=''
        self.bandData=[]
        self.dossData=[]
        self.allBands=[]
        self.bandx=[]
        self.dpi=150
        self.title=''.join(re.split('(\.)',title)[:-2])
        self.outFileName=title
        self.bandFileName=self.title+'_dat.BAND'
        self.dossFileName=self.title+'_dat.DOSS'
        self.bandRange='10000,10000'
        self.combined='1'
        self.mass='n'
        self.intersect=0
        self.gaussian=None
        self.pdosColumn=[]
        self.pdosLabel=[]
        self.lineColors=[]
        self.pathData=[]
        self.labelk=''
        self.kIndex=''
        self.KLabelList=[]
        self.KLabelListNum=[]
        self.KLabelIndex=[]
    def readOut(self):
        outFile=open(self.outFileName)
        for line in outFile:
            if re.search('\*\s+FROM BAND\s+\d+\s+TO BAND',line):
                self.startBand=int(line.split()[3])
                self.lastBand=int(line.split()[6])
            if re.search('N. OF ELECTRONS PER CELL',line):
                self.HOMO=int(line.split()[5])/2
            if re.search('FERMI ENERGY\s+-?\d+',line):
                self.fermi=float(line.split()[line.split().index('FERMI')+2])
        outFile.close()

    
    def readBand(self):
        bandRawData=open(self.bandFileName)
        for line in bandRawData:
            if re.search('-?\d+E[-+]\d+',line):
                self.bandData.append(line.split())
            if re.search('^#\s+\d+\s+',line):
                self.KLabelIndex.append(line.split()[1])
            if re.search('^@ XAXIS TICKLABEL\s+\d+',line):
                self.KLabelList.append(line.split()[3].split('"')[1])
            if re.search('^@ XAXIS TICK\s+\d+',line):
                #print line.split()
                self.KLabelListNum.append(line.split()[4])
        #read k point label and index
        #for kpoint in self.pathData:
        #    kplist=map(float,filter(None,re.split("\D+",kpoint[1])))
        #    kplabel=''
        #    print kplist
        #    for i in range(3):
        #        if kplist[i]/kplist[3] == 0:
        #            kplabel=kplabel+'0-'
        #        else:
        #            kplabel=kplabel+'{0:.2f}-'.format(kplist[i]/kplist[3])[1:]
        #    self.KLabelList.append(kplabel[:-1])
        #    self.KLabelIndex.append(kpoint[0])
        bd=pylab.array(self.bandData)
        bandDATA=bd.transpose()
        #cut a segment of k path
        if self.kIndex:
            firstK=int(self.kIndex.split(',')[0])
            secondK=int(self.kIndex.split(',')[1])
            newBandData=[]
            for i in bandDATA:
                newBandData.append(i[self.KLabelIndex[firstK]:int(self.KLabelIndex[secondK])+1])
            bandDATA=newBandData
        
        
        #sort bandDATA to get rid of crossing
        sortedbd=[]
        self.bandX=map(float,bandDATA[0])
        bd=pylab.transpose(bandDATA[1:])
        for line in bd:
            fline=map(float,line)
            sortedbd.append(list(sorted(fline)))
        bandDATA=[list(i+self.fermi) for i in pylab.transpose(sortedbd)]
        #print bandDATA
        self.bandDATA=bandDATA
        self.bandXmin=float(self.bandX[0])
        self.bandXmax=float(self.bandX[-1])
        self.bandYmin=min(bandDATA[0])
        self.bandYmax=max(bandDATA[-1])
        self.xpoints=len(bandDATA[0])

        #self.bandYmin=pylab.array(map(float,bd.transpose()[1])).min()+self.fermi
        #self.bandYmax=pylab.array(map(float,bd.transpose()[-1])).max()+self.fermi
        #check for band accross and restore order
        if self.intersect==1:
            '''check for band accross and restore order
            first find points where two bands are close to each other (contact points)
            then search the the contact point to find the switch point that can make the slope of the bands most smooth
            finally sort the bands according to max of occ. band and min of virt. band
            
            '''
            for i in range(len(bandDATA)-2):
                for j in range(2,len(bandDATA[i])-2):     
                    for k in range(i+1,len(bandDATA)-1):
                        if pylab.absolute(bandDATA[i][j]-bandDATA[k][j]) < 0.1:
                            #while pylab.absolute(bandDATA[i][j]-bandDATA[k][j]) > pylab.absolute(bandDATA[i][j+1]-bandDATA[k][j+1]):
                            #    j=j+1
                            if pylab.absolute(2*bandDATA[i][j]-bandDATA[i][j-1]-bandDATA[i][j+1]) > pylab.absolute(2*bandDATA[k][j]-bandDATA[i][j-1]-bandDATA[k][j+1]):
                                while pylab.absolute(2*bandDATA[k][j]-bandDATA[i][j-1]-bandDATA[k][j+1]) > pylab.absolute(2*bandDATA[k][j+1]-bandDATA[i][j]-bandDATA[k][j+2]) and j+3 < len(bandDATA[k]):
                                    j=j+1
                                topl=bandDATA[k][:j]
                                topr=bandDATA[k][j:]
                                btnl=bandDATA[i][:j]
                                btnr=bandDATA[i][j:]
                                bandDATA[i]=btnl+topr
                                bandDATA[k]=topl+btnr
            sortedOcc=sorted(bandDATA[:self.HOMO],key=lambda kpoint: max(kpoint))
            sortedVirt=sorted(bandDATA[self.HOMO:],key=lambda kpoint: min(kpoint))
            bandDATA=sortedOcc+sortedVirt  
        self.allBands=bandDATA
        
         
    def readDOSS(self):
        try: dossRawData=open(self.dossFileName)
        except IOError:
            dossRawData=''
            return
        for line in dossRawData:
            if re.search('-?\d+E[-+]\d+',line):
                self.dossData.append(map(float,line.split()))
    
        sortedDoss=sorted(self.dossData,key=lambda l:l[0])
        self.dossData=sortedDoss
        self.npdos=len(self.dossData[0])-1
        self.dossX=pylab.array(map(float,pylab.array(self.dossData).transpose()[0] ))
        self.pdossY=[]
        try:
            if max(map(int,self.pdosColumn)) > self.npdos-1:
                print 'number of input pdos column '+str(self.pdosColumn)+' exceed number of pdos in file: '+str(self.npdos-1)
        except ValueError:
            pass
        for i in range(1,self.npdos):
            if str(i) in self.pdosColumn:
                self.pdossY.append(pylab.array(map(float,pylab.array(self.dossData).transpose()[i] )))
        self.dossY=pylab.array(map(float,pylab.array(self.dossData).transpose()[-1] ))
        if self.gaussian>0:
            self.FWHM=((self.dossX[0]-self.dossX[-1])/len(self.dossX)*self.gaussian)*27.211396*2.355
            WdossY=gaussian_filter1d(self.dossY,self.gaussian,truncate=100.0,mode='constant')
            WpdossY=[]
            for i in range(len(self.pdossY)):
                WpdossY.append(gaussian_filter1d(self.pdossY[i],self.gaussian,truncate=100.0,mode='constant'))
            self.dossY=WdossY
            self.pdossY=WpdossY
        self.dossYmin=self.dossY.min()
        self.dossYmax=self.dossY.max()
        
    def propGen(self):
        #calculate bandwidth and bandGap for HOMO and LUMO 
        HOMO=self.HOMO-self.startBand
        self.HOMOband=pylab.array(self.allBands[HOMO])
        self.LUMOband=pylab.array(self.allBands[HOMO+1])
        HOMOmax=self.HOMOband.max()
        HOMOmin=self.HOMOband.min()
        LUMOmax=self.LUMOband.max()
        LUMOmin=self.LUMOband.min()
        self.HOMOmax=HOMOmax*27.211396
        self.LUMOmin=LUMOmin*27.211396
        self.HOMObandWidth=(HOMOmax-HOMOmin)*27.211396
        self.LUMObandWidth=(LUMOmax-LUMOmin)*27.211396
        self.bandGap=(LUMOmin-HOMOmax)*27.211396
        if 'h' in self.mass.lower():
            #calculate effective mass
            '''1 cut the band into n segment
            the segment involve the HOMO max or LUMO min will be used for fitting
            use parabole to fit the curve, calculate effective mass, R square and mean absolute error
            if the R squre > 0.9995 and MAE < 0.0001. the effective mass will be adopted'''
            fitErrorHOMO=[]
            fitRHOMO=[]
            massHole=[]
            allCoefficients=[]
            partHOMOs=[]
            partXHs=[]
            fitHOMOs=[]
            fitX=[]
            for n in range(2,self.xpoints/5):
                maxIndex=self.HOMOband.argmax()
                if maxIndex<len(self.HOMOband)/n:
                    partHOMO=self.HOMOband[:len(self.HOMOband)/n]
                    partXH=self.bandX[:len(self.HOMOband)/n]
                elif maxIndex>len(self.HOMOband)/n*(n-1):
                    partHOMO=self.HOMOband[-len(self.HOMOband)/n:]
                    partXH=self.bandX[-len(self.HOMOband)/n:]
                else:
                    partHOMO=self.HOMOband[maxIndex-len(self.HOMOband)/(2*n):maxIndex+len(self.HOMOband)/(2*n)]
                    partXH=self.bandX[maxIndex-len(self.HOMOband)/(2*n):maxIndex+len(self.HOMOband)/(2*n)]

                coefficients=pylab.polyfit(partXH,partHOMO,2)
                polynomial = pylab.poly1d(coefficients)
                d2=polynomial.deriv().deriv()
                massHole.append(d2(self.bandX[maxIndex])**-1)
                fitHOMO=polynomial(partXH)
                HOMObar=pylab.sum(partHOMO)/len(partHOMO)
                ssreg=pylab.sum((fitHOMO-HOMObar)**2)
                sstot=pylab.sum((partHOMO-HOMObar)**2)
                fitErrorHOMO.append(pylab.sum(pylab.absolute(fitHOMO-partHOMO))/len(partHOMO))
                fitRHOMO.append(ssreg/sstot)
                fitX.append(n) 
                #following four list for drawing the fit graph
                allCoefficients.append(coefficients)
                partHOMOs.append(partHOMO)
                partXHs.append(partXH)
                fitHOMOs.append(fitHOMO)
            self.holeMass='---'
            minError=1
            minErrorIndex=1
         
            for n in range(self.xpoints/5-2):
                if fitRHOMO[n]>0.9995:
                    if minError > fitErrorHOMO[n]:
                        minError=fitErrorHOMO[n]
                        minErrorIndex=n
            self.holeMass=massHole[minErrorIndex] 
            fitX=[i**-1 for i in fitX]
            pylab.figure('homofit')
            pylab.plot(partXHs[minErrorIndex],partHOMOs[minErrorIndex],'r-')
            pylab.plot(partXHs[minErrorIndex],fitHOMOs[minErrorIndex],'b-',label=r'$E={0[0]:.3G}k^2{0[1]:+.3G}k{0[2]:+.3G}$'.format(allCoefficients[minErrorIndex])+ \
            '\n'+r'$m=(\frac{\partial^2E}{\partial k^2})^{-1}=$'+r'${:.3G} m_0$'.format(self.holeMass)+
            '\n'+r'$R^2={:.4G}$'.format(fitRHOMO[minErrorIndex]) + \
            '\n'+r'$MAE={:.4G}$'.format(fitErrorHOMO[minErrorIndex]) +\
            '\n$fraction={:.2G}$'.format((minErrorIndex+2)**-1))
            pylab.legend()
            pylab.ylabel(r'Energy (hartree)')
            pylab.xlabel(r'$k  (bohr^{-1})$')
            pylab.savefig(self.title+'_HoleMassFit.svg',dpi=self.dpi,bbox_inches='tight',format='svg')
            pylab.clf()
            pylab.figure(3)
            pylab.subplot(3,1,1)
            pylab.plot(fitX,massHole,'b-')
            pylab.ylabel('Effective Mass')
            pylab.subplot(3,1,2)
            pylab.plot(fitX,fitRHOMO,'b-')
            pylab.ylabel('R squre')
            pylab.subplot(3,1,3)
            pylab.plot(fitX,fitErrorHOMO,'b-')
            pylab.ylabel('Mean Absolute Error')
            pylab.xlabel('fraction of the k path')
            pylab.savefig(self.title+'_massFit.png',dpi=self.dpi,bbox_inches='tight',format='svg')
        else:
            self.holeMass='---'
        if 'e' in self.mass.lower():
            #calculate effective mass
            '''1 cut the band into n segment
            the segment involve the HOMO max or LUMO min will be used for fitting
            use parabole to fit the curve, calculate effective mass, R square and mean absolute error
            if the R squre > 0.9995 and MAE < 0.0001. the effective mass will be adopted'''
            fitErrorLUMO=[]
            fitRLUMO=[]
            massElectron=[]
            allCoefficients=[]
            partLUMOs=[]
            partXLs=[]
            fitLUMOs=[]
            fitX=[]
            
            for n in range(2,self.xpoints/5):
                minIndex=self.LUMOband.argmin()
                if minIndex<len(self.LUMOband)/n:
                    partLUMO=self.LUMOband[:len(self.LUMOband)/n]
                    partXL=self.bandX[:len(self.LUMOband)/n]
                elif minIndex>len(self.LUMOband)/n*(n-1):
                    partLUMO=self.LUMOband[-len(self.LUMOband)/n:]
                    partXL=self.bandX[-len(self.LUMOband)/n:]
                else:
                    partLUMO=self.LUMOband[minIndex-len(self.LUMOband)/(2*n):minIndex+len(self.LUMOband)/(2*n)]
                    partXL=self.bandX[minIndex-len(self.LUMOband)/(2*n):minIndex+len(self.LUMOband)/(2*n)]
                coefficients=pylab.polyfit(partXL,partLUMO,2)
                polynomial = pylab.poly1d(coefficients)
                d2=polynomial.deriv().deriv()
                fitLUMO=polynomial(partXL)
                LUMObar=pylab.sum(partLUMO)/len(partLUMO)
                ssreg=pylab.sum((fitLUMO-LUMObar)**2)
                sstot=pylab.sum((partLUMO-LUMObar)**2)
                fitErrorLUMO.append(pylab.sum(pylab.absolute(fitLUMO-partLUMO))/len(partLUMO))
                fitRLUMO.append(ssreg/sstot)
                massElectron.append(d2(self.bandX[minIndex])**-1)
                fitX.append(n)  
                #following four list for drawing the fit graph
                allCoefficients.append(coefficients)
                partLUMOs.append(partLUMO)
                partXLs.append(partXL)
                fitLUMOs.append(fitLUMO)
            self.electronMass='---'
            minError=100
            minErrorIndex=0
            for n in range(self.xpoints/5-2):
                if fitRLUMO[n]>0.999:
                    if minError==0:
                        minError=fitErrorLUMO[n]
                    elif minError > fitErrorLUMO[n]:
                        minError=fitErrorLUMO[n]
                        minErrorIndex=n
            self.electronMass=massElectron[minErrorIndex] 
            fitX=[i**-1 for i in fitX]
            pylab.figure('lumofit')
            pylab.plot(partXLs[minErrorIndex],partLUMOs[minErrorIndex],'r-')
            pylab.plot(partXLs[minErrorIndex],fitLUMOs[minErrorIndex],'b-',label=r'$E={0[0]:.3G}k^2{0[1]:+.3G}k{0[2]:+.3G}$'.format(allCoefficients[minErrorIndex])+ \
            '\n'+r'$m=(\frac{\partial^2E}{\partial k^2})^{-1}=$'+r'${:.3G} m_0$'.format(self.electronMass)+
            '\n'+r'$R^2={:.4G}$'.format(fitRLUMO[minErrorIndex]) + \
            '\n'+r'$MAE={:.4G}$'.format(fitErrorLUMO[minErrorIndex])+\
            '\nr$fraction={:.2G}$'.format(minErrorIndex**-1)
            )
            pylab.legend()
            pylab.ylabel(r'Energy (hartree)')
            pylab.xlabel(r'$k  (bohr^{-1})$')
            pylab.savefig(self.title+'_ElectronMassFit.png',dpi=self.dpi,bbox_inches='tight',format='png')
            pylab.figure(3)
            pylab.subplot(3,1,1)
            pylab.plot(fitX,massElectron,'r-')
            pylab.ylabel('Effective Mass')
            pylab.subplot(3,1,2)
            pylab.plot(fitX,fitRLUMO,'r-')
            pylab.ylabel('R squre')
            pylab.subplot(3,1,3)
            pylab.plot(fitX,fitErrorLUMO,'r-')
            pylab.ylabel('Mean Absolute Error')
            pylab.savefig(self.title+'_massFit.png',dpi=self.dpi,bbox_inches='tight',format='png')
            pylab.xlabel('fraction of the k path')
        else:
            self.electronMass='---'
        
        
    def plotBand(self):
        downRange=int(self.bandRange.split(',')[0])
        upRange=int(self.bandRange.split(',')[1])
        HOMO=self.HOMO-self.startBand
        startBand=HOMO-downRange
        lastBand=HOMO+upRange
        if startBand < 0:
            startBand=0
        if lastBand > self.lastBand-self.startBand:
            lastBand=self.lastBand-self.startBand
        occBand=self.HOMO-self.startBand
        virBand=lastBand-HOMO
        pylab.figure(1)
        if not self.dossData:
            self.combined='0'
        if self.combined=='1':
            pylab.subplot(121)
        for i in range(startBand,lastBand+1):
            if i <= occBand:
                pylab.plot(self.bandX,pylab.array(self.allBands[i])*27.211396,'b-')
            if i > occBand:
                pylab.plot(self.bandX,pylab.array(self.allBands[i])*27.211396,'r-')  
        
        if self.labelk == '1':
            pylab.xticks([self.bandX[int(i)-1] for i in self.KLabelIndex],self.KLabelList)
        elif self.labelk == 'n':
            pylab.xticks([self.bandX[int(i)-1] for i in self.KLabelIndex],self.KLabelListNum)
        elif self.labelk:
            customLabelList=self.labelk.upper().split(',')
            if len(customLabelList) != len(self.KLabelIndex):
                print " the number of inputed K label "+str(len(customLabelList))+" is not equal to number of K points in file:"+str(len(self.KLabelIndex))
                sys.exit()
            pylab.xticks([self.bandX[int(i)-1] for i in self.KLabelIndex],customLabelList)
             
        bandYmin=pylab.array(self.allBands[startBand]).min()*27.211396
        bandYmax=pylab.array(self.allBands[lastBand]).max()*27.211396
        bandYmargin=(bandYmax-bandYmin)*0.05
        pylab.axis([self.bandXmin,self.bandXmax,bandYmin-bandYmargin,bandYmax+bandYmargin])
        pylab.ylabel('Energy (eV)')
        pylab.xlabel('k (Bohr-1)')
        pylab.minorticks_on()
        if self.combined=='0':
            pylab.savefig(self.title+'_bands.png',dpi=self.dpi,bbox_inches='tight',format='png')
            pylab.savefig(self.title+'_bands.svg',bbox_inches='tight',format='svg')
        if self.combined=='1':
            pylab.subplot(122)
        if self.combined=='0' and self.dossData:
            pylab.figure(2)
        #plot dos and pdos 
        if self.dossData:
            pylab.plot(self.dossY,self.dossX*27.211396,'b-')
            if self.pdosColumn:
                
                for i in range(len(self.pdossY)):
                    pylab.plot(self.pdossY[i],self.dossX*27.211396,color=self.lineColors[i],label=self.pdosLabel[i])
                pylab.legend()
            pylab.axis([self.dossYmin,self.dossYmax,bandYmin-bandYmargin,bandYmax+bandYmargin])
            #write doss date to file
            dossDATA=open('DOSSDATA','w')
            for i in range(0,len(self.dossX)):
                dossDATA.write('%10f %f\n' % (self.dossX[i],self.dossY[i]*4))
            dossDATA.close()
        if self.combined=='0' and self.dossData:
            pylab.ylabel('Energy (eV)')
        pylab.xlabel('State Density')
        pylab.minorticks_on()
        if self.combined=='1':
            pylab.savefig(self.title+'_band-doss.png',dpi=self.dpi,bbox_inches='tight',format='png')
            pylab.savefig(self.title+'_band-doss.svg',bbox_inches='tight',format='svg')
        if self.combined=='0' and self.dossData:
            pylab.savefig(self.title+'_doss.png',dpi=self.dpi,bbox_inches='tight',format='png')
            pylab.savefig(self.title+'_doss.svg',bbox_inches='tight',format='svg')
        pylab.figure(2)
        pylab.clf()
        pylab.figure(1)
        pylab.clf()
parser=argparse.ArgumentParser(description='Generate band structure figure and corresponding data')
parser.add_argument("-b","--band",metavar='bandFile',help='Band structure file generate by crystal14, default:outfile_dat.BAND',type=str)
parser.add_argument("-d","--doss",metavar='dossFile',help='Dstructure file generate by crystal14, default:outfile_dat.DOSS',type=str)
parser.add_argument("-r","--range",help='range for draw band structure, e.g. -1,2 means draw form HOMO-1 to LUMO+1, default:all the band',type=str)
parser.add_argument("-m","--mass",help='fit the band structure to calculate effective mass,h=hole, e=electron,other=not calculate, default:0',default='n',type=str)
parser.add_argument("-c","--combined",help='draw Band structure and Doss seperately or combined,0=combined, 1=seperately, default:1',type=str)
parser.add_argument("-s","--size",help='the dpi of the figure,0=combined, 1=seperately, default:1',type=str)
parser.add_argument("-g","--gaussian",help='gaussian widen for the DOS, set the sigma, default:None',type=int )
parser.add_argument("-p","--pdos",help='input three tuple. The first one is the column to plot, the second one is the Label, and the third one is the style, e.g. 1,2,3;S,Px,Py;r,g,b, default:None',type=str )
parser.add_argument("-l","--labelk",help='input labels for k points; seperated by comma, e.g. g,k,m. input 1 for fraction k points:e.g. 0.5,0,0', type=str)
parser.add_argument("-k","--kIndex",help='set the k point index to draw band structure, and to fit effective mass, e.g. 0,1 use first and second k points', type=str)
parser.add_argument("outputFile",nargs='+')
args=parser.parse_args()
print "%-35s%-12s%-12s%-12s%-15s%-15s%-12s%-12s" % ('file name','HOMO Width','LUMO Width','Band Gap','Hole Mass','Electron Mass','HOMO','LUMO')

def genColorGrad(step):
    rg=step/2+1
    gb=step/2+1
    if rg+gb>step+1:
        gb=gb-1
    rglist=[]
    gblist=[]

    for i in range(rg):
        rglist.append((1-i/(rg-1.0),i/(rg-1.0),0.0))
    for i in range(gb):
        gblist.append((0.0,1.0-i/(rg-1.0),i/(rg-1.0)))
    rglist.pop()
    return rglist+gblist

for outfile in args.outputFile:
    bandgen=bandProp(outfile)
    if args.band:
        bandgen.bandFileName=args.band
    if args.labelk:
        bandgen.labelk=args.labelk
    if args.doss:
        bandgen.dossFileName=args.doss
    if args.range:
        bandgen.bandRange=args.range
    if args.combined:
        bandgen.combined=args.combined
    if args.kIndex:
            bandgen.kIndex=args.kIndex
    if args.pdos:
        print args.pdos
        bandgen.pdosColumn=args.pdos.split(';')[0].split(',')
        try:
            bandgen.pdosLabel=args.pdos.split(';')[1].split(',')
        except IndexError:
            bandgen.pdosLabel=['']*len(bandgen.pdosColumn)
        if len(bandgen.pdosLabel)< len(bandgen.pdosColumn):
            bandgen.pdosLabel=bandgen.pdosLabel+[' ']*(len(bandgen.pdosColumn)-len(bandgen.pdosLabel))
        try:
            bandgen.lineColors=args.pdos.split(';')[2].split(',')
        except IndexError:
            bandgen.lineColors=genColorGrad(len(bandgen.pdosColumn))
        #print bandgen.lineColors
    bandgen.gaussian=args.gaussian
    bandgen.mass=args.mass
    bandgen.readOut()
    bandgen.readBand()
    bandgen.readDOSS()
    bandgen.propGen()
    print "%-35s%-12.6f%-12.6f%-12.6f%-15.6s%-15.6s%-12.6f%-12.6f" % (bandgen.title,bandgen.HOMObandWidth,bandgen.LUMObandWidth,bandgen.bandGap,bandgen.holeMass,bandgen.electronMass,bandgen.HOMOmax,bandgen.LUMOmin)
    bandgen.plotBand()
    if bandgen.gaussian>0:
        print bandgen.title+': FWHM is '+str(bandgen.FWHM)+'(eV)'
    #print bandgen.HOMOband
    
    
    #if last_max==pylab.array(map(float,line)).min()
    #last_max=pylab.array(map(float,line)).max()
    #last_min=pylab.array(map(float,line)).min()
#    





